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ABSTRACT 

We study the dynamical stability against bar-mode deformation of rapidly spinning neutron stars with 
differential rotation. We perform fully relativistic 3D simulations of compact stars with M/R > 0.1, 
where M is the total gravitational mass and R the equatorial circumferential radius. We adopt an 
adiabatic equation of state with adiabatic index r = 2. As in Newtonian theory, we find that stars above 
a critical value of (3 = T/W (where T is the rotational kinetic energy and W the gravitational binding 
energy) are dynamically unstable to bar formation. For our adopted choices of stellar compaction and 
rotation profile, the critical value of (3 = (3dGR is ~ 0.24 — 0.25, only slightly smaller than the well-known 
Newtonian value ~ 0.27 for incompressible Maclaurin spheroids. The critical value depends only very 
weakly on the degree of differential rotation for the moderate range we surveyed. All unstable stars form 
bars on a dynamical timescale. Models with sufficiently large (3 subsequently form spiral arms and eject 
mass, driving the remnant to a dynamically stable state. Models with moderately large (3 > f3d.GR do not 
develop spiral arms or eject mass but adjust to form dynamically stable ellipsoidal-like configurations. 
If the bar-mode instability is triggered in supernovae collapse or binary neutron star mergers, it could 
be a strong and observable source of gravitational waves. We determine characteristic wave amplitudes 
and frequencies. 



1. INTRODUCTION 

Neutron stars in nature are rotating and subject to non- 
axisymmetric rotational instabilities. An exact treatment 
of these instabilities exists only for incompressible equilib- 
rium fluids in Newtonian gravity (see, e.g., Chandrasekhar 
1969; Tassoul 1978; Shapiro & Teukolsky 1983). For these 
configurations, global rotational instabilities arise from 
nonradial toroidal modes e ml(p (m = ±1,±2,...) when 
(3 = T/W exceeds a certain critical value. Here ip is the 
azimuthal coordinate and T and W are the rotational ki- 
netic and gravitational potential binding energies. In the 
following we will focus on the m = ±2 bar mode, since it is 
the fastest growing mode when the rotation is sufficiently 
rapid. 

There exist two different mechanisms and correspond- 
ing timescales for bar mode instabilities. Uniformly ro- 
tating, incompressible stars in Newtonian theory are secu- 
larly unstable to bar mode formation when (3 > (3 S ~ 0.14. 
However, this instability can only grow in the presence of 
some dissipative mechanism, like viscosity or gravitational 
radiation, and the growth time is determined by the dis- 
sipative timescale, which is usually much longer than the 
dynamical timescale of the system. By contrast, a dy- 
namical instability to bar mode formation sets in when 
(3 > (3d — 0.27. This instability is independent of any dis- 
sipative mechanisms, and the growth time is determined 
by the hydrodynamical timescale of the system. 

The secular instability in compressible stars, both uni- 
formly and differentially rotating, has been analyzed nu- 
merically within linear perturbation theory, by means of 
a variational principle and trial functions, by solving the 
eigenvalue problem, or by other approximate means. This 



technique has been applied not only in Newtonian the- 
ory (Lynden-Bell & Ostriker 1967; Ostriker & Boden- 
heimer 1973; Ipser and Lindblom 1989; Friedman and 
Schutz 1977) but also in post-Newtonian theory (Cutler 
and Lindblom 1992; Shapiro and Zane 1998 for incom- 
pressible stars) and full general relativity (Yoshida and 
Eriguchi 1995; Bonazzola, Frieben and Gourgoulhon 1996; 
Stergioulas and Friedman 1998). For relativistic stars, the 
critical value of [3 S depends on the compaction M/R of the 
star (where M is the gravitational mass and R the circum- 
ferential radius at the equator), on the rotation law and 
on the dissipative mechanism. The gravitational-radiation 
driven instability occurs for smaller rotation rates, i.e. for 
values (3 S < 0.14, in general relativity. For extremely com- 
pact stars (Stergioulas and Friedman 1998) or strongly dif- 
ferentially rotating stars (Imamura et al. 1995), the criti- 
cal value can be as small as (3 S < 0.1. By contrast, viscos- 
ity drives the instability to higher rotation rates (3 S > 0.14 
as the configurations become more compact (Bonazzola, 
Frieben and Gourgoulhon 1996; Shapiro and Zane 1998). 

Determining the onset of the dynamical bar-mode in- 
stability, as well as the subsequent evolution of an unsta- 
ble star, requires a numerical simulation of the fully non- 
linear hydrodynamical equations. Simulations performed 
in Newtonian theory (e.g. Tohline, Durisen & McCol- 
lough 1985; Durisen, Gingold & Tohline 1986; Williams 
& Tohline 1988; Houser, Centrella & Smith 1994; Smith, 
Houser & Centrella 1995; Houser & Centrella 1996; Pick- 
ett, Durisen & Davis 1996; New, Centrella & Tohline 1999) 
have shown that (3d depends only very weakly on the stiff- 
ness of the equation of state. Once a bar has developed, 
the formation of spiral arms plays an important role in 
redistributing the angular momentum and forming a core- 
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halo structure. Recently, it has been shown that, similar 
to the onset of secular instability, fid can be smaller for 
stars with a higher degree of differential rotation (Tohline 
& Hachisu 1990; Pickett, Durisen & Davis 1996) 

To date, the dynamical bar-mode instability has been 
analyzed only in Newtonian theory, hence almost nothing 
is known about the role of relativistic gravitation. The 
reason is that until quite recently a stable numerical code 
capable of performing reliable hydrodynamic simulations 
in three dimensions plus time in full general relativity has 
not existed. Some recent developments, however, have 
advanced the field significantly. New formulations of the 
Einstein equation based on modifications of the standard 
3 + 1 ADM system of equations (Arnowitt, Deser & Mis- 
ner 1962) have resulted in codes which have proven to be 
remarkably stable over many dynamical timescales (e.g., 
Shibata & Nakamura 1995; Baumgarte & Shapiro 1999; 
Oohara & Nakamura 1999). In addition, gauge conditions 
which warrant long-time stable evolution for rotating and 
self-gravitating systems and are manageable computation- 
ally have been developed (e.g., Shibata 1999b). In this 
paper we adopt the relativistic hydrodynamic implemen- 
tation of Shibata (1999a) to study the onset and growth 
of the dynamical bar-mode instability in relativistic stars. 
Although this study is carried out only for a simple equa- 
tion of state and rotational law, it demonstrates how, as 
numerical relativity in full 3 + 1 matures, it is becoming 
more useful as a tool to solve long-standing problems in 
relativistic astrophysics characterized by strong gravita- 
tional fields and little or no spatial symmetry. 

There are numerous evolutionary paths which may lead 
to the formation of rapidly rotating neutron stars with 
[3 ~ 0.3. The parameter (3 increases approximately as 
during stellar collapse. During supernova collapse, 
the core contracts from ~ 1000 km to ~ 10 km, and hence 
p increases by about two orders of magnitude. Thus, even 
moderately rapidly rotating progenitor stars may yield 
rapidly rotating neutron stars which may reach the on- 
set of dynamical instability (Bonazzola & Marck 1993; 
Rampp, Muller & Ruffert 1998). Similar arguments hold 
for accretion induced collapse of white dwarfs to neutron 
stars and for the merger of binary white dwarfs to neu- 
tron stars. In fact, recent X-ray and radio observations 
of supernova remnants have identified several young, iso- 
lated, rapidly rotating pulsars, suggesting that these stars 
may have been born with periods of several milliseconds 
(Marshall et al. 1998; Kaspi et al. 1998; Torii et at. 1999). 
These neutron stars could be the collapsed remnants of 
rapidly rotating progenitors. 

Rapidly rotating neutron stars may naturally arise in 
the merger of binary neutron stars. Baumgarte, Shapiro 
and Shibata (2000) have studied equilibrium configura- 
tions of differentially rotating neutron stars and found ex- 
amples where the maximum allowed mass increases by a 
factor of about 2 due to differential rotation. This sug- 
gests that the merger of binary neutron stars could result 
in a "hypermassive" neutron star which has rest mass ex- 
ceeding the maximum value for uniformly rotating stars. 
Recent hydrodynamic simulations in full general relativity 
indicate that such hypermassive neutron stars can indeed 
be produced in the merger of moderately compact neu- 
tron stars (Shibata & Uryu 2000). They show that the 
remnant is unlikely to exceed the onset point of dynami- 



cal instability initially. Subsequent neutrino emission and 
cooling, however, will make the star shrink in size, leading 
to an increase in /3, possibly beyond the onset of nonradial 
dynamical instability, (3d- 

Rapidly rotating neutron stars experiencing the bar- 
mode instability could have significant observable conse- 
quences. According to Newtonian simulations (Tohline, 
Durisen & McCollough 1985; Durisen Gingold & Tohline 
1986; Williams & Tohline 1988; Houser, Centrella & Smith 
1994; Smith, Houser & Centrella 1995; Houser & Centrella 
1996), a dynamically unstable star may evolve into a two- 
component system containing a central star and circum- 
stellar accretion disk. Such a system may be observable 
in a supernova remnant. In the case of merged binaries, 
the differentially rotating remnant may be more massive, 
hot and bloated than a typical rapidly rotating, old pul- 
sar. Consequently, the frequency of gravitational waves 
excited by the bar- mode instability could be low, i.e., less 
than 1kHz (see Eq. ( pc| ) below), and hence detectable by 
kilomctcr-sizc laser interferometers such as LIGO (Lai & 
Shapiro 1995; Thorne 1995). 

In this paper we summarize the results of our fully rel- 
ativistic simulations of bar-mode instabilities in neutron 
stars. We determine (3d for highly relativistic stars, follow 
the growth of the bar-mode, and find the frequency and 
amplitude of the emitted gravitational waves. We imple- 
ment the numerical scheme described in Shibata (1999a), 
using differentially rotating neutron stars of high (3 for ini- 
tial data. We focus on differentially rotating stars, since 
uniformly rotating stars do not reach (3 ;> 0.2 except for 
extremely stiff equations of state, and hence do not be- 
come dynamically unstable to bar- modes (Tassoul 1978). 
We adopt an adiabatic equation of state with T = 2 as 
a reasonable qualitative approximation to moderately stiff 
nuclear equation of state. The adiabatic assumption is jus- 
tified even for hot neutron stars, since energy dissipation 
is small over the dynamical timescales of interest. 

In Sec. 2, we briefly summarize our formulation of the 
fully relativistic system of equations and our numerical 
scheme. In Sec. 3, initial models of differentially rotating, 
equilibrium neutron stars are presented. Following Shi- 
bata, Baumgarte & Shapiro (2000) we adopt the so-called 
conformal flatness approximation to prepare differentially 
rotating neutron stars in (approximate) equilibrium states 
for computational convenience. To confirm the reliability 
of this approximation, we also compute numerically exact 
equilibrium states and demonstrate that this approxima- 
tion is accurate (cf. Cook, Shapiro & Teukolsky 1996). 
In Sec. 4, we present our numerical results, focusing on 
the onset of the bar-mode instability, its early growth and 
corresponding waveforms and frequencies. We briefly sum- 
marize our results in Sec. 5. 

Throughout this paper, we adopt geometrized units with 
G = 1 = c where G and c denote the gravitational constant 
and speed of light. In numerical simulation, we use Carte- 
sian coordinates x k = (x, y, z) with r = \J x 2 + y 2 + z 2 , 
w = \/ x 2 + y 2 and ip = tan _1 (y/x); t denotes coordinate 
time. Greek indices /i, v, . . . denote x, y, z and t, and Latin 
indices k, . . . denote x, y and z. 



2. SUMMARY OF THE FORMULATION 
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We perform hydro-dynamic simulations in full 3 + 1 gen- 
eral relativity (GR). We use the same formulation and 
gauge conditions as in Shibata (1999a), to which the reader 
may refer for details and basic equations. The fundamen- 
tal variables used in this paper are: 

p : rest mass density, 
e : specific internal energy, 
P : pressure, 
: four velocity, 



a 
I3 k 
7ij 

7 



lapse function, 
shift vector, 

metric in 3D spatial hypersurface, 
= e 12 ^ = det( 7ij ), 



Kn : extrinsic curvature. 



Geometric variables, (f>, 7^, the trace of the extrinsic cur- 
vature K = K^, Aij = e'^iKij -HjK/3), as well 
as three auxiliary functions F-i = djjij, where dj is the 
partial derivative, are evolved with an unconstrained evo- 
lution code in a modified form of the ADM formalism (Shi- 
bata & Nakamura 1995). GR hydrodynamic equations are 
evolved using a van Leer scheme for the advection terms 
(van Leer 1977; Hawley, Smarr & Wilson 1984). Numer- 
ical simulation is performed using Cartesian coordinates. 
Violations of the Hamiltonian constraint and conservation 
of mass and angular momentum are monitored as code 
checks. Several test calculations, including spherical col- 
lapse of dust, stability of spherical neutron stars, and the 
stable evolutions of rigidly and rapidly rotating neutron 
stars have been described in Shibata (1999a). Simulations 
using this code and exploring the dynamical (quasi-radial) 
stability against gravitational collapse of rigidly rotating 
"supramassive" neutron stars, which have rest masses ex- 
ceeding the TOV limit for a nonrotating spherical star, 
have been presented in Shibata, Baumgarte & Shapiro 
(2000). A simulation using this code and demonstrating 
the existence of dynamically stable differentially rotating 
"hypermassive" stars, which have rest masses exceeding 
the maximum value for uniformly rotating stars, was pre- 
sented in Baumgarte, Shapiro and Shibata (2000). 
The stress energy tensor for an ideal fluid is given by 



(p + pe + P)u,j,u v + Pg^ 



(1) 



where g^ v is the spacetime metric. We adopt a T-law 
equation of state 



P = (T - l)pe, 



(2) 



where T is the adiabatic constant. For isentropic configu- 
rations the T-law equation of state can be rewritten in the 
polytropic form 



P = np L 



1 

n 



(3) 



where k is the polytropic constant and n the polytropic 
index. This is the form that we use for constructing initial 



data. Throughout this paper, we adopt n — 1 as a rea- 
sonable qualitative approximation to a moderately stiff, 
nuclear equation of state for simplicity. 

Instead of p and e we numerically evolve the densities 
= pau°e 6 ^ and e* = (pe) 1 / r o;u e 6 ^ as the hydrody- 
namic variables (Shibata, Oohara & Nakamura 1997; Shi- 
bata 1999a). Since these variables satisfy evolution equa- 
tions in conservation form, the total rest mass of the sys- 
tem 



M = 



xp* 



(4) 



is automatically conserved, as is the the volume integral 
of the energy density e* in the absence of shocks. 

The time slicing and spatial gauge conditions we use 
in this paper for the lapse and shift are the same as 
those adopted in our series of papers (Shibata 1999a, 
1999b; Shibata, Baumgarte & Shapiro 2000); i.e. we im- 
pose an "approximate" maximal slice condition (K ~ 0) 
and an "approximate" minimum distortion gauge condi- 
tion (Di(d t j l: >) — where Di is the covariant derivative 
with respect to 7^ , see Shibata 1999b). 

3. INITIAL CONDITIONS FOR ROTATING NEUTRON STARS 

As initial conditions, we adopt rapidly and differen- 
tially rotating neutron stars in (approximate) equilibrium 
states. The approximate equilibrium states are obtained 
by choosing a conformally flat spatial metric, i.e., assum- 
ing 7y = e^Sij (see, e.g., Cook, Shapiro & Teukolsky, 
1996, or Shibata 1999a for the equations to be solved in 
this approximate framework). This approach is computa- 
tionally convenient and, as demonstrated in Cook, Shapiro 
& Teukolsky (1996) and Shibata, Baumgarte & Shapiro 
(2000) provides an excellent approximation to exact ax- 
isymmetric equilibrium configurations in rigid rotation. 

Following previous studies (e.g. Komatsu, Eriguchi & 
Hachisu 1989a, 1989b; Cook, Shapiro & Teukolsky 1992, 
1994; Bonazzola, Gourgoulhon, Salgado & Marck 1993; 
Salgado, Bonazzola, Gourgoulhon, & Haensel 1994; Gous- 
sard, Haensel & Zdunik 1998) we fix the differentially ro- 
tational profile according to 



= u°u v = A 2 (n - O), 



(5) 



where A is an arbitrary constant with dimensions of length 
(which describes the length scale over which £1 changes) 
and Qo is the angular velocity on the rotation axis, which 
is chosen to be the z-axis. In the Newtonian limit u° — > 1 
and u v — > n7 2 0, the rotational profile reduces to 



A 2 n Q 

+ A 2 ' 



(6) 



Thus, for smaller A, Q becomes a steeper function of w. 

Equilibrium configurations of rotating stars are charac- 
terized by their gravitational mass M, angular momentum 
J, rotational kinetic energy T and gravitational potential 
binding energy W, which in GR can be defined invariantly 
according to 



M 



= /<- 



2T n ° + T^)ae 6 *d 3 x, 



(7) 
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J = / T*ae*+&x 



(8) 

T = r / rJT, o °ae 60 d 3 a;, (9) 



>3 



PF = / p*edTx + T + M - M. 



(10) 



(e.g. Cook, Shapiro & Teukolsky 1992). For approximate 
configurations derived in the conformal flatness approxi- 
mation, on the other hand, the gravitational mass is com- 
puted from the asymptotic behavior of the conformal fac- 
tor, which, after using Gauss' law and the Hamiltonian 
constraint yields 



M 



1 



(p + pe + P){av°) 2 - P + —K ij K,, 



16tt 



o e 



5 *d 3 x 



(11) 

(see, e.g., Bowen & York 1980). This expression is correct 
independent of axisymmetry. For conformally flat configu- 
rations, J, T, and W can be computed from Eqs. (8)-(10) 
for nonaxisymmetric configurations as well. As in Newto- 
nian gravity, we define (3 as the ratio T/W. (Note W > 
in our definition.) 

Physical units enter the problem only through the poly- 
tropic constant k, which can be chosen arbitrarily or else 
completely scaled out of the problem. In the following, we 
define 

M = M aK - n/2 , M = Mn- n/2 , 



J 

Pmax 



J K , P„ 



Prot K 



-n/2 



(12) 



where P rot and p max are rotational period and the maxi- 
mum density. Note that p max does not necessarily coincide 
with the central density for stars of highly differential rota- 
tion. The barred quantities are now independent of k, and 
all results can be scaled for arbitrary k using Eqs. (|l2|). 

For the construction of (approximate) equilibrium mod- 
els for initial data, we adopt a grid in which the semi-major 
axes of the stars, chosen along the x and y-axes in the equa- 
torial plane, are covered with 40 grid points. The semi- 
minor axis in the polar direction along z-axis is covered 
with ~ 10 grid points for the case (3 ~ 0.25 and A ~ r e . 
Hereafter, the coordinate lengths of the semi-major and 
minor axes are referred as r e and r p , respectively. We 
have confirmed the convergence of our numerical solutions 
by increasing the number of grid points covering r e to 120 
for typical models shown in Table I. Comparing with these 
higher resolution models shows that the numerical error in 
M, T, W and J of the lower resolution models is less than 
1 %. 

We find that for n = 1, stars with (3 ;> 0.2 can only 
be constructed for A <; r e . Also, for A < r e /2, stars 
with j3 ;> 0.25 are not spheroid but toroids. Focusing 
on spheroidal stars with f3 ^> 0.2, we study cases with 

A = A/r e = 0.8 and 1. 

In Table I, we display parameters for selected models 
constructed both within the conformal flatness approxi- 
mation and from the exact numerical equations (using the 

3 In the conformal flatness approximation, the error in T and W is 
1996). Hence, the magnitude of the error for T and W is expected to 
~ (GM / Rc 2 )^ < 1%, which is consistent with our numerical results. 



code from Cook, Shapiro & Teukolsky 1994) for a given 
Pmax, A and r p /r e . Note that both approaches lead to 
valid, fully relativistic initial data in the sense that both 
data satisfy the constraint equations of Einstein's field 
equations. The difference is that the "exact" solutions 
provide an exactly stationary solution, while the confor- 
mally flat solutions may only be approximately station- 
ary when evolved dynamically. Here, we choose stars of 
R/M ~ 7 and 9.5 (Hereafter, R denotes the circumfer- 
ential radius of the equator). All the stars in Table I 
are "hypermassive" (Baumgarte, Shapiro & Shibata 2000) 
with rest masses larger than the maximum rest mass of 
rigidly rotating stars built from the same equation of state 
(M ~ 0.207; cf. Fig. 2). The matter profiles of the equi- 
librium configurations obtained in the conformal flatness 
approximation agree fairly well with the exact numerical 
solution. The slight deviation arises mainly from the error 
associated with the conformal flatness approximation (as 
opposed to errors due to the different finite differencing 
in the two different codes). As can be seen in Table I, 
P is systematically underestimated in the conformal flat- 
ness approximation by < 0.004 (i.e., < 2%) although the 
deviation for Mq and M is considerably less (< 1% er- 
ror). From a post-Newtonian point of view, T and W are 
quantities of 0(c~ 2 ) but Mo and M are of O(c ). In the 
conformal flatness approximation, we neglect the second 
post-Newtonian terms of 0(c~ 4 ) in the metric (cf. Kley & 
Schafer 1999) so that the error for T and W can be larger 
than that for masses by 0(c 2 ) 3 . Indeed, the agreement 
between exact and approximate solutions is improved for 
less compact stars (compare models D2 and D3 with D6 
and D7). Because of this small deviation from the exact 
solution, the initial conditions prepared in the conformal 
flatness approximation should be regarded as slightly per- 
turbed states of exact equilibria. 

In Fig. 1, we plot Q/Qq as a function of w in the equa- 
torial plane to illustrate the rotational velocity field for 
different A. We show an unstable (to bar-modes) star 
with R/M - 7 and (3 ~ 0.25, for both A = 1 and 0.8. 
As expected from Eq. (^|), Q is a steeper function of w 
for smaller A. As demonstrated in the Newtonian calcula- 
tions of Pickett, Durisen and Davis (1996), the degree of 
falloff of £1 versus w could be an important factor for de- 
termining the onset point of nonaxisymmetric dynamical 
instability in differentially rotating stars. 

4. NUMERICAL RESULTS 

To investigate the dynamical stability against bar-mode 
deformation, we initially superimpose a density perturba- 
tion of the form 



P = Po ( 1 + $b 



y 



(13) 



where po denotes the axisymmetric configuration, and we 
choose 5b to be 0.1 or 0.3. We leave the four- velocity Ui 
unperturbed. We recompute the constraint (initial value) 
equations whenever we modify the equilibrium configura- 
tions this way, to guarantee that we are satisfying the Ein- 
stein equations at t = 0. 

0(c~ 4 ), and 0(c~ e ) in M and Mq (see, e.g., Asada and Shibata, 
be ~ (GM/Rc 2 ) 2 ~ 2% for R/M = 7, but that for M and M is 
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The growth of a bar mode can be followed by monitoring 
the distortion parameter 



V 



where x z lms denotes the mean square axial length 



Mo J 



d 3 xp* (x l f 



1/2 



(14) 



(15) 



For dynamically unstable stars, r\ grows exponentially un- 
til reaching a saturation point, while for stable stars, it 
remains approximately constant for many rotational peri- 
ods. 

We perform simulations using a fixed uniform grid with 
typical size 153 x 77 x 77 in x — y — z and assume 7r-rotation 
symmetry around the z-axis as well as a reflection symme- 
try about the z = plane. We have also performed test 
simulations with different grid resolutions to check that 
the results do not change significantly. Since we impose 
7r-rotation symmetry, we ignore one-armed spiral (m = 1) 
modes which might be dominant for rotating stars in which 
is a very steep function of w (Pickett, Durisen & Davis 
1996). However, 7r-rotation symmetry guarantees that the 
center-of-mass drift is identically zero (cf. New, Centrella 
& Tohline 1999). 

We note that the outer boundaries of our computational 
domain reside inside the wavelength of gravitational waves 
emitted by the bar- mode perturbation, A gw . The typ- 
ical location of the outer boundaries along each axis is 
~ (0.1 — 0.2)A gw . Without setting the boundaries in the 
radiation zone at r ;> A gw , or using a sophisticated wave- 
extraction technique in the near zone (Bishop et al. 1996; 
Abrahams et al. 1998), it is impossible to compute the ra- 
diation reaction and asymptotic waveforms completely ac- 
curately. However, in this paper we focus on dynamical in- 
stabilities, which are independent of dissipation processes 
and grow on a dynamical timescale considerably shorter 
than the secular dissipation timescale due to gravitational 
wave emission. The error in the evaluation of the gravita- 
tional waves can therefore be safely neglected in assessing 
the onset and growth of a dynamical instability. 

Fig. 2 summarizes our findings on the dynamical stabil- 
ity against bar-mode deformation. We evolve a range of 
stellar models with A = 0.8 (squares) and A = 1 (circles). 
All of these models have (3 > 0.2 and 0.2 < r p /r e < 0.35, 
and we determine their stability by inducing an initial non- 
axisymmetric perturbation 5b = 0.1. Each model takes 
50 - 100 CPU hours to run on the FACOM VX/4R ma- 
chine; the long-time runs described below take about 150 
CPU hours. Stable stars are denoted with open circles and 
squares, and unstable stars with solid circles or squares. 
For models denoted by crosses we were unable to deter- 
mine stability unambiguously. Note that all differentially 
rotating stars shown here are dynamically stable against 
quasi-radial collapse to black holes. 

Fig. 2 shows that in a Mo versus p ma x diagram an un- 
stable region can be clearly separated from a stable region. 
The demarcation line is nearly independent of the degree 
of differential rotation, at least for the modest variation in 
the rotation law that we consider (recall that for (3 ;> 0.2 

spheroidal stars only exist within a restricted range of A) . 



The two regions are separated by a thick dashed-dotted 
line in the figure. We have also plotted lines of constant 
(3; one where (3 = 0.245 for A = 1 and another where 
(3 = 0.24 for A = 0.8. These two lines very closely trace 
the demarcation line between the regions of stability and 
instability. This result suggests that in general relativity 
as in Newtonian gravitation, the parameter (3 is a good 
diagnostic for assessing whether a rotating star is stable 
against the dynamical bar mode instability. It also suggests 
that for differentially rotating, relativistic stellar models, 
the threshold for dynamical bar formation (3dGR may de- 
pend only weakly on the differential rotation law, and is 
only slightly smaller than the corresponding value for uni- 
formly rotating, Newtonian stars, (3d ~ 0.27. 

The small but measurable decrease in the critical value 
of f3 could be either due to the presence of differential ro- 
tation (cf. Tohline & Hachisu (1990) and Pickett, Durisen 
& Davis (1996), who have observed this effect in Newto- 
nian gravity), or by GR effects (compare with the decrease 
in (3 S with increasing compaction for the secular onset of 
the gravitational-wave driven instability in Stcrgioulas & 
Friedman 1998), or a combination of both. In order to 
separate the two effects it would be desirable to system- 
atically explore parameter space and study models with 
different rotation laws and varying compaction up to the 
Newtonian limit of small M/R. We are preparing such a 
survey now (Saijo, Shibata, Baumgarte & Shapiro 2000). 

In anticipation of this survey, recall that in Newto- 
nian gravity, Mo scales with pmax 1 ^ 2 ™ (Shapiro & Teukol- 
sky 1983) for polytropes, or Mo oc p max for n = 1. 
For stars with f3 ~ 0.27, this relation turns out to be 
M ~ 10 — 12p max ; for smaller values of (3 the coefficient 
is only slightly smaller. We therefore expect that the line 
marking the onset of dynamical instability, a line of con- 
stant f3, approaches a linear relationship in Fig. 2 in the 
Newtonian limit. 

Unfortunately, using a fully relativistic code is im- 
practical for simulating stars in the Newtonian or post- 
Newtonian regime. The Courant condition restricts the 
numerical timestep to the light travel time across a 
grid zone, and therefore scales with R. The dynamical 
timescale of the star, however, is approximately the free- 
fall timescale i? 3 / 2 /M 1 / 2 , so that the ratio between the dy- 
namical timescale and Courant timestep scales roughly as 
[R/M) 1 / 2 . In the Newtonian limit this ratio becomes very 
large, so that many timesteps have to be carried out in or- 
der to simulate a fixed number of dynamical timescales, 
which makes the calculation computationally impracti- 
cal. To avoid this problem, we arc implementing a post- 
Newtonian code (Shibata, Baumgarte & Shapiro 1998) to 
explore the intermediary regime, and will present these re- 
sults in a forthcoming paper (Saijo, Shibata, Baumgarte 
& Shapiro 2000). Our preliminary finding is that in the 
Newtonian limit, (3d ~ 0.26 for A = 1, and hence the onset 
of dynamical instability occurs slightly earlier for stars of 
greater compaction. 

For five of the above models (Dl, D2, D3, D6 and D7) 
we have followed the growth and evolution of the bar-mode 
instability over several rotational timescales. In Figs. 3, 4 
and 5 we show snapshots of density contours and velocity 
fields in the x — y (left) and x — z (right) planes for the 
compact models Dl, D2 and D3 {R/M ~ 7; see Table I). 
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The models are rotating counterclockwise. Note that the 
stars start out as highly flattened, disk-like objects. In 
Fig. 6, we also show r\ as a function of time (a) for models 
Dl, D2 and D3 as well as (b) for models D6 and D7, which 
are slightly less compact (R/M ~ 9.5). In order to accel- 
erate the growth of the instabilities, we set 5b = 0.3 for 
these simulations. In addition, we plot the early evolution 
of \r]\ on a semi-log scale in Fig. 7. This plot also demon- 
strates that while the perturbation parameter 5b — 0.3 
is fairly large, the nondimensional measure |r/| is initially 
safely in the linear regime. Unstable growth ceases once 
the bar-mode reaches nonlinear saturation. Following the 
evolution much beyond this point is impractical, both be- 
cause of accumulation of numerical error, and because fur- 
ther evolution begins to be affected by gravitational wave 
emission, which is crudely handled in this code as discussed 
above. 

Model Dl is stable against bar- mode formation and 
demonstrates the ability of our code to identify and main- 
tain such a configuration in stable equilibrium (see Fig. 3). 
The other four models are unstable and the barlike per- 
turbation (and hence rj) grows exponentially in the early 
phase until saturation is reached, typically when n ~ 
0.2—0.4. However, the evolution after the saturation varies 
for different models. For model D2, spiral arms form in 
the outer part of the barlike object, which then spread 
outward, transporting away mass and angular momentum 
(see Fig. 4). In Fig. 8, we show the fraction of the rest 
mass inside a fixed coordinate radius r, M*(r)/M*, as a 
function of time for model D2. We find that a few per cent 
of the total rest mass is ejected from the star. We also find 
that the fraction of the rest mass inside r = 0.2r e and 0.6r e 
ultimately increases during the late phase of the evolution 
since the mass with high specific angular momentum is 
transported outwards, and the star slightly contracts. As 
a consequence, the maximum density of the star increases 
at late times by about 50% from p max — 0.045 to ~ 0.065 
by t = 6P r " t 4 . Here P® ot is the rotational period on the ro- 
tation axis (w = 0). Tracking its motion in Fig. 2, we find 
that the unstable star approaches the stability threshold 
line from the left, and ultimately enters the stable region. 
Model D6 evolves very similarly. These findings are also 
in qualitative agreement with Newtonian results (Tohline, 
Durisen & McCollough 1985; Durisen Gingold & Tohline 
1986; Williams & Tohline 1988; Houser, Centrella & Smith 
1994; Smith, Houser & Centrella 1995; Houser & Centrella 
1996; New, Centrella & Tohline 1999). 

Unstable models D3 (see Fig. 5) and D7 start out closer 
to the stability threshold ((3 ;> PdGFt) and never form spi- 
ral arms or eject mass. Instead, they evolve to an ellip- 
soidal shape, which is maintained for many rotational pe- 
riods. Although the mass is not ejected in this case, the 
maximum density slightly increases again due to outward 
angular momentum transport. Consequently, the stars 
again approach the threshold line of the dynamical sta- 
bility shown in Fig. 2 and become stable. 

We determine the growth time and oscillation period of 
the bar-mode instability from Figs. 7 and summarize the 
results in Table II. From Figs. 7, we can match r\ to a 



function 

T] ~ TfalO* 77 " 9 cos(27rf/T + ip ) (16) 

where r/ and <f are constants, and r g and r Q are the 
growth time and oscillation period. The growth time r g 
depends strongly on the compaction R/M and (3, and is 
smaller for larger (3 as expected. The oscillation period 
t ~ 1.2P r " t depends only very weakly on the compaction 
R/M and (3, and approximately agrees with that found 
in Newtonian simulations for n — 1 (Williams & Tohline 
1988). The characteristic oscillation period of n after the 
saturation of the growth is r Q <~ (1.2 — 1.4)P r ° t for all 
models, so that the pattern period is <~ (2.4 — 2.8)P? ot . 

In order to check the numerical convergence of our re- 
sults, we repeated these simulations with a lower resolu- 
tion (101 x 51 x 51 as opposed to 153 x 77 x 77, with the 
outer boundaries at the same location). In Figs. 6 (a) and 
(b), the dotted lines denote the low-resolution result for 
rj as a function of time. For t/P? ot <^ 5 the results agree 
well, implying that a fair qualitative convergence has been 
achieved. For later times, the accumulation of numerical 
truncation error and problems associated with the inade- 
quate outer boundaries results in a poorer agreement be- 
tween the two resolutions. 

The differentially rotating ellipsoids formed after sat- 
uration are highly flattened and still seem to have high 
(3 ;> 0.2. We expect these to be secularly unstable against 
gravitational radiation (Stergioulas & Friedman 1998). 
Therefore, rj will probably maintain a fairly high value 
of O(0. 1) on a radiation reaction timescale, which is much 
longer than the rotational period. As argued by Lai & 
Shapiro (1995), such an ellipsoid will ultimately settle 
down to an axisymmetric star or a Dcdckind-likc ellipsoid 
(a nonaxisymmetric, stationary star whose figure does not 
rotate but which has internal differential motion, cf. Chan- 
drasekhar 1969). 

As a measure of gravitational waveforms we show h + 
and h x 5 

h+ = r{% x - % y )/2M (17) 
h x = r% y /M (18) 

in Fig. 9 for models Dl (dotted lines), D2 (solid lines) and 
D3 (dashed lines). These quantities are read off near the 
outer boundaries along the z-axis as a function of retarded 
time. For models D6 and D7 very similar waveforms to 
those of D2 and D3 are generated, as expected from Fig. 6. 
As mentioned above, the outer boundaries are not located 
in the wave zone, which implies that these quantities will 
not exactly agree with the asymptotic waveforms. We 
guess that the wave shapes are in fairly good agreement 
with the exact ones because the frequency agrees with the 
oscillation frequency of the bar pattern, but the error in 
the amplitude may be ;> 10% or larger. 

We find that for unstable stars, the maximum magni- 
tude of h + and h x is ~ 0.03 — 0.08 depending on (3 and 
R/M, while they remain of O(10 -3 ) for the stable star 
Dl. The maximum observable wave amplitude from such 
a source situated at a distance r from the earth is then 



4 For the stable model Dl we find a small increases in the central density as well, but only by about 10 % after t = 6P r " t . This small increase 
is due to numerical viscosity and is a numerical artifact. The central density in the model D2 increases by a much larger fraction over the same 
timescale, suggesting that the bulk of this increase is indeed a physical effect. 

5 Our quantities /i+ and h x differ from those defined in Misner, Thorne & Wheeler (1973) by a factor r/M. 
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approximately 



ft, ~ 6 x 10" 



-22 / "+,X 



0.05 / V 2.5M 



The typical gravitational wavelength after satura- 
tion is ~ 1.3P r Q t or, using the empirical relation 
P r Q ot - 3M(R/M) 3 / 2 (cf. Table I), approximately ~ 
70M(i?/7M) 3/2 . This implies that the frequency of grav- 
itational waves is 



V 7M / V M J 



(20) 



We note that the amplitude and frequency are in ap- 
proximate agreement with earlier Newtonian calculations 
(Houser, Centrella & Smith 1994; Smith, Houser & Cen- 
trella 1995; Houser & Centrella 1996), which suggests that 
GR effects do not drastically alter the simple quadrupole- 
formula predictions for the waveforms. 

Eqs. ( |l9| ) and (20) give the maximum amplitude for one 
cycle ana its frequency. As discussed in Lai & Shapiro 
(1995), the effective amplitude can be much larger because 
of the quasi-periodic nature of the source. Also, the fre- 
quency will gradually shift to smaller values as a result of 
radiation reaction. Thus, even if h in one cycle is small 
and the frequency is initially as high > 1kHz, these nearly- 
ellipsoidal stars may eventually be observable by kilome- 
ter size laser-interferometric gravitational wave detectors 
like LIGO (Thorne 1995) for stars with M ~ 2.5M Q and 
R ~ 7M as a result of the secular decrease of the fre- 
quency. Clearly, this predicted drift needs to be confirmed 
by a more detailed study. Also, gravitational waves could 
be a source for a specially designed narrow band interfer- 
ometers or resonant-mass detectors in which the frequency 
band is between 1 and 2 kHz (Thorne 1995). 

5. SUMMARY 

We have performed numerical simulations of rapidly and 
differentially rotating neutron stars in full 3+1 general rel- 
ativity. We treated compact stars of 10 p> R/M ^6 and 
focussed on their dynamical stability against bar-mode for- 
mation. We found that when plotted in a Mq versus p max 
diagram, a region of stable stars can be clearly distin- 
guished from a region of unstable stars, with the onset 
of instability almost independent of the degree of differen- 
tial rotation. We showed that the parameter (3 = T/W 
remains a good diagnostic of the onset point of insta- 
bility in the relativistic domain as it did for Newtonian 
stars. The critical value for the instability onset depends 



only weakly on the degree of differential rotation for the 
models surveyed to date. For those cases we find that 
PdGR ~ 0.24 — 0.25, and that (3dGR decreases slightly for 
stars with a higher degree of differential rotation. We 
also have preliminary evidence that f3dGR decreases with 
compaction as well. We will systematically study this hy- 
pothesis with a post-Newtonian numerical analysis in a 
forthcoming paper (Saijo, Shibata, Baumgarte & Shapiro 
2000). 

For selected models, we followed the growth and satu- 
ration of bar-mode perturbations up to late times. Stars 
with sufficiently large (3 > f3d.GR develop bars first and 
then form spiral arms, leading to mass ejection. Stars 
with smaller values of (3 ~ /3dGR also develop bars, but do 
not form spiral arms and eject only very little mass. In 
both cases, unstable stars appear to form differentially ro- 
tating, triaxial ellipsoids once the bar-mode perturbation 
saturates. Typically, these flattened ellipsoids appear to 
have (3 f> 0.2, so that they would be secularly unstable 
due to gravitational waves and viscosity. We expect that 
this secular instability will allow the stars to maintain a 
bar-like shape for many dynamical timescales, leading to 
quasi-periodic emission of gravitational waves. 

We estimate the initial frequency and amplitude of grav- 
itational waves to be / ~ (1 — 1.4)kHz and h ~ 5xl0 -22 for 
stars of mass ~ 2.5M and radius R ~ 7M at a distance 
of lOMpc. The effective amplitude of gravitational waves 
could be much larger during the subsequent evolution be- 
cause of the accumulation of quasi-periodic wave cycles 
(Lai & Shapiro 1995). In order to accurately determine 
the secular evolution of the ellipsoidal star together with 
emitted gravitational wave signal, a more detailed calcu- 
lation is necessary. Since the secular timescale is larger 
than the dynamical timescale by several orders of magni- 
tude, it may be impossible to follow the evolution with 
a fully dynamical code, even with implicit differencing to 
avoid the Courant criterion for stability. This suggests 
that in full GR, the secular evolution problem may best 
be solved within an appropriate, quasi-stationary scheme 
similar in spirit to the approach used in stellar evolution 
calculations. 

Numerical computations were performed on the VX/4R 
machines in the data processing center of the National 
Astronomical Observatory of Japan. This work was sup- 
ported by NSF Grants AST 96-18524 and PHY 99-02833 
and NASA Grant NAG5-7152 at the University of Illinois 
at Urbana-Champaign (UIUC). M.S. gratefully acknowl- 
edges support by JSPS (Fellowships for Research Abroad) 
and the hospitality of the Department of Physics at UIUC. 
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Table I. Maximum density p max , rest mass Mo, gravitational mass M, compaction R/M, (3 = T/W, rotation period on 
the rotation axis w = (-P r " t ) and at the equator (Pf ot ) of differentially rotating stars for selected models computed in 
the conformally flatness approximation (upper line) and in the exact equations (lower line). For Dl - D5, R/M ~ 7 and 
for D6 and D7, it is ~ 9.5. J/M 2 is larger than unity for the models shown here. 



1 pi 1 e 


A 


n 

pmax 


iv±o 


M 


R/M 


T/W 


pa 
-'rot 


-'rot 


stahilitv 


Model 


0.35 


1 


06056 


0.260 


0.241 


6.62 


0.230 


12.7 


36.7 


tD u < J > s I V 


Dl 








0.259 


0.241 


6.61 


0.233 


12.7 








0.275 


1 


0.04460 


0.277 


0.259 


7.07 


0.258 


15.0 


41.8 


unstable 


D2 








0.277 


0.258 


7.08 


0.262 


14.9 








0.30 


1 


0.04590 


0.264 


0.246 


7.24 


0.251 


14.9 


41.3 


unstable 


D3 








0.264 


0.246 


7.23 


0.254 


14.8 








0.25 


0.8 


0.04650 


0.262 


0.245 


7.02 


0.243 


12.3 


45.3 


unstable 


D4 








0.261 


0.244 


7.02 


0.247 


12.2 








0.325 


0.8 


0.05940 


0.254 


0.235 


6.50 


0.223 


10.5 


40.2 


stable 


D5 








0.253 


0.235 


6.50 


0.226 


10.4 








0.275 


1 


0.03070 


0.229 


0.217 


9.27 


0.262 


19.8 


50.7 


unstable 


D6 








0.229 


0.217 


9.27 


0.265 


19.8 








0.3 


1 


0.03220 


0.219 


0.208 


9.41 


0.254 


19.5 


49.7 


unstable 


D7 








0.219 


0.208 


9.39 


0.256 


19.4 









Table II. Growth time T g and oscillation period r Q of the bar-mode for unstable stars D2, D3, D6, and D7. The timescales 
are shown in units of P" t . 



Model 


T a 


To 


D2 


2.48 


1.21 


D3 


3.62 


1.19 


D6 


2.07 


1.23 


D7 


2.85 


1.21 
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Fig. 1. — n/fio as a function of zu in the equatorial plane for differentially rotating stars of R/M ~ 7 and (3 ~ PdGR and for A = 1, 0.8 
and 2/3. 
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Fig. 2. — Models of differentially rotating stars in a Mo versus pmax diagram. Circles denote stars with A — 1, scjuarcs with A — 0.8. Solid 
(open) circles or squares represent stars that are unstable (stable). Marginally stable stars are denoted with a cross. The region for the stable 
stars is clearly separated from that for the unstable stars by the thick dashed-dotted line. This line is followed fairly closely by the dashed and 
dotted lines, which have been constructed for differentially rotating stars of (A, /3) = (1,0.245) and (0.8,0.24). We carry out long-duration 
simulations for those stars denoted by a hyphen (models Dl, D2, D3, D6 and D7). The long-dashed and solid lines are for nonrotating 
spherical stars and rigidly rotating stars at the mass shedding limit. Scales for the top and right axes are shown for re = 100(G 3 Mg/c 4 ) in 
which the maximum rest mass for spherical stars is about 1.8Mq. 
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Fig. 3. — Snapshots of density contours for p» and the velocity flow for v l in the equatorial plane (left) and in the y = plane (right) for 
the stable model Dl. The contour lines are drawn for p*/p* max = 10~ 0,3j for j = 0, 1, 2, • ■ • , 10 where p* max is 0.193, 0.209 and 0.248 for 
the three different times (the corresponding values of p max are 0.061, 0.064, and 0.067). The lengths of arrows are normalized to 0.3c (left) 
and 0.1c (right). The time is shown in units of P? ot - 



SHIBATA, BAUMGARTE & SHAPIRO 



13 



t = 0.00P 




-10 -5 5 10 2 4 6 8 10 

X / M X / M 



t = 5.82P 




-10 -5 5 10 2 4 6 8 10 

X / M X / M 



Fig. 4. — Snapshots of density contours for p» and the velocity flow for v % in the equatorial plane (left) and in the y = plane (right) for 
the unstable model D2. The contour lines are drawn for p*/p» max = 10 _0 - 3j for j = 0, 1, 2, • • ■ , 10 where p* max is 0.126, 0.172 and 0.264 for 
the three different times (the corresponding values of p m ax are 0.045, 0.059, and 0.065). The lengths of arrows are normalized to 0.3c (left) 
and 0.1c (right). The time is shown in units of P? ot - 
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FlG. 5. — Snapshots of density contours for p* and the velocity flow for v % in the equatorial plane (left) and in the y = plane (right) for 
the unstable model D3. The contour lines are drawn for p* / p* max — 

l0-0-3j f or j = 0, 1, 2, • • • , 10 where p* max is 0.128, 0.152 and 0.176 for 
the three different times. The lengths of arrows are normalized to 0.3c (left) and 0.1c (right). The time is shown in units of P^ ot - 
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Fig. 6. — The distortion parameter tjasa function of t (a) for models Dl (dashed line), D2 (solid line) and D3 (long-dashed line), and (b) 
for models D6 (solid line) and D7 (long-dashed line). The results in a low resolution simulation with 101 X 51 X 51 grid points are shown by 
the dotted lines. 




Fig. 7. — \r)\ as a function of t (a) for models D2 (solid line) and D3 (dotted line), and (b) for models D6 (solid line) and D7 (dotted line). 
The dottcd-dashed lines denote the growth time of the bar-mode instability. 
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Fig. 9. — h+ and h x as a function of a retarded time t — z ^ s for stars Dl (dotted line), D2 (solid line) and D3 (dashed line). 



